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ABSTRACT 

This paper aims at quantifying discreetness effects, born of finite particle number, on 
the dynamics of dark matter haloes forming in the context of cosmological simulations. 
By generalising the standard calculation of two body relaxation to the case when the 
size and mass distribution are variable, and parametrising the time evolution using 
established empirical relations, we find that the dynamics of a million particle halo is 
noise-dominated within the inner percent of the final virial radius. Far larger particle 
numbers (~ 10 s ) are required for the RMS perturbations to the velocity to drop to 
the 10% level there. The radial scaling of the relaxation time is simple and strong: 
^roiax ~ f j implying that numbers 3> 10 s are required to faithfully model the very in- 
ner regions; artificial relaxation may thus constitute an important factor, contributing 
to the contradictory claims concerning the persistence of a power law density cusp to 
the very centre. The cores of substructure haloes can be many relaxation times old. 
Since relaxation first causes their expansion before recontraction occurs, it may render 
them either more difficult or easier to disrupt, depending on their orbital parameters. 
It may thus modify the characteristics of the subhalo distribution; and, if as suggested 
by several authors, it is parent-satellite interactions that determine halo profiles, the 
overall structure of the system may be affected. We derive simple closed form formulas 
for the characteristic relaxation time of both parents and satellites, and an elementary 
argument deducing the weak iV-scaling reported by Diemand et al. (2004) when the 
main contribution comes from relaxing subhaloes. 

Key words: dark matter - galaxies: haloes - diffusion - galaxies: general - galaxies: 
formation - galaxies: structure 



1 INTRODUCTION 

-/V-body modelling of the dynamical evolution of the dark 
matter component involves the approximation whereby the 
system to be simulated is represented by a surrogate one 
consisting of many order of magnitudes fewer particles. Dis- 
creteness noise is thus necessarilly greatly enhanced. Fur- 
thermore, within the context of the cold dark matter sce- 
nario, material rapidly collapses into haloes; in a hierarchi- 
cal process that ensures that the first structures are badly 
resolved, with particle noise propagating through the hier- 
archy, further enhancing discreteness effects and weakening 
the iV-dependence of the effective relaxation time (Binney 
& Knebe 2002; Diemand et al. 2004). 

Relaxation is expected to play its most prominent role 
in the very central regions of collapsed structures — pre- 
cisely the place where there remains considerable contro- 
versy as to whether the slope of the density profile persists as 
a power law up to the resolution limit (Diemand et al. 2005) , 
or instead gradually flattens into a smooth core (Navarro 



et al. 2004). From the difference between these two situa- 
tions follow important implications concerning the compat- 
ibility of the cold dark matter scenario with the observed 
inner mass distribution in galaxies. Relaxation also affects 
the symmetry of a system in physical and velocity space, 
rendering it more isotropic; and triaxial haloes may play an 
important role in determining the dynamics of the baryonic 
component (e.g., El-Zant & Shlosman 2002; El-Zant et al. 
2003). 

Since substructure haloes are resolved with far fewer 
particles, relaxation is expected to have a more dominant 
effect on their internal structure and (by modifying the 
efficiency of stripping) their spatial distribution. If accre- 
tion and merging are are dominant processes in determining 
them (Syer & White 1998; Dekel et al. 2003b; Dekel et al. 
2003a), this will have important implications for the parent 
halo properties too. 

Estimating the importance of particle noise in iV-body 
simulations is an essential step in evaluating the role it may 
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play in determining the dynamical state of collapsed struc- 
tures. Our focus here is on the generalisation of the standard 
formulation of the problem, leading to the familiar formulas 
for the 'relaxation time' in stellar systems, to the cosmolog- 
ically relevant case, when the physical extent and the mass 
distribution of the object under consideration continuously 
vary with time. The details of the procedure, along with the 
assumptions involved, are outlined in Section[5] In our treat- 
ment we will distinguish between 'parent' haloes, represent- 
ing the most massive progenitors and 'subhaloes', represent- 
ing the substructure. The former are assumed to continue to 
grow in mass and size, while the latter's growth stops once 
they are incorporated into a larger structure. In both cases 
we derive simple formulas for calculating characteristic re- 
laxation times (Section|3J and, in addition, undertake a more 
detailed calculation of the expected root mean square (RMS) 
perturbation to trajectories' velocities for haloes evolving in 
the context of the 'concordance' CDM model; giving sim- 
ple empirical fits to the results (Sectional . Conclusions are 
summarised in Section 



2 THE MODEL 

We are interested in the description of relaxation due to 
discreteness (particle) noise in evolving cosmological haloes. 
The model through which this is done is described below. 
First the time evolution of the halo mass distribution is 
parametrised; the standard description of relaxation in grav- 
itational systems is then generalised to encompass the case 
of these evolving systems. 



2.1 Evolution of the halo mass distribution 

Describing the temporal evolution of the halo mass distribu- 
tion is simplified by the fact that their densities can always 
be approximately parametrized by double power law func- 
tions, independently of the mass of the halo, or the redshift 
at which it is identified. Two profiles have been extensively 
used: that introduced by Navarro, Frenk & White 1997 
(NFW) has density ~ 1/r in the inner region and ~ 1/r 3 in 
the outer one; the associated mass fraction within radius r 



M(r) _ ln(l + r/r s ) 



r jr s 
1 + r jr s 



M(r vir ) 



m(l + c)-^- 



(1) 



The profile put forward by Moore et al. 1999 (M99) has 
similar asymptotic density variation at large r, but a steeper 
central cusp (~ l/r 3/2 ). It has mass distribution 



M(r) _ ln[l + (r/r a ) 3/2 ] 



M(r vir ) ln(l + c 3 / 2 ) 



(2) 



In these formulas, r s is the characteristic halo scale length 
(separating the 1/r from the 1/r 3 region) and r v i r is the 
virial radius. Their ratio c = r v i r /r s determines the mean 

density contrast inside radius r: A — jjjy-^-ff- (cf. Fig0. 

The virial radius is usually defined as the radius inside 
which the average density is ~ 200 times the critical density 



for closure p c . Relative to p c , the mean density contrast 
inside r is 

„=£M = 200 xA = 200 -^l% (3) 

pc M(r vir ) r 3 v ' 

In a universe with matter and vacuum energy densities 
such that Q m (z) + Qy(z) = 1, the virial radius varies with 
redshift as 



»WO) _ (M(r vir ,z)/M ) 1/3 



(4) 



where M = M(r vir ,z = 0). Wechsler et al. (2002) sug- 
gested the empirical relation whereby 

M vil (z) = M e- az . (5) 
This defines a characteristic halo growth time 1 



/dlnMV 
V dt ) 



adz I dt 



(6) 



The parameter a ~ 1 varies systematically with halo mass 
(but appears constrained to the range 0.4 < a < 1.6). In this 
paper we fix a — 1, noting that our results were found to 
be quite insensitive as to the exact choice of a. Wechsler et 
al. also find that, for NFW type fits to the density profile, 
c ~ (1 + with an average value of order 10 at z = 0. 

Their results thus generally confirmed the prescription (due 
to Bullock et al. 2001) 



[NFW, z] = 



(7) 



which is used here. We assign a minimal value of c = 9/5 
at high z (because structures with p ~ 1/r in their outer 
regions are necessarily far from virial equilibrium; the simu- 
lations of Bullock et al. 2001; Zhao et al. 2003 and Tasitsiomi 
et al. 2004 confirm the presence of a minimal c). The con- 
centration parameter of an M99 fit relates to to that of an 
NFW fit through c[M99] w (c[NFW]/1.7) 9 . 

For any given z there is a mass dependent spread in c; 
it is given, in terms of the typical mass scale M* (z) collaps- 
ing at redshift z, via c ~ (M v i r /M,) _0 13 . Nevertheless, this 
relation shows that, even for a halo that is ten orders of mag- 
nitude more massive than M* , c only changes by a factor of 
ten. From Fig. 1, this leads to a corresponding change in the 
density contrast of about the same factor or less. But since 
the relaxation time (cf, Eq. ll5l below'l ~ N(r)rD(r) ~ \J p(r) 
inside any fixed radius; it changes by a factor of only a few 
even for such a highly unlikely mass deviation. Trials with 
different normalisations in Eq. 0, and with the prescription 
proposed by (Eq. 10 of) Zhao et al. (2003), have confirmed 
this. In fact it turns out (Section 4.2) that the relaxation rate 
is quite a strong function of radius, and will be important 
at small radii for any reasonable set of c values during the 
evolution. The results presented in this paper assume the 



1 The exponential form reflects a two phased evolution; a period 
of rapid accretion followed by slower growth. Zhao et al. (2003) 
have confirmed its basic premises for halos growing through al- 
most four orders of magnitudes in mass and three in radius (al- 
though they do not use this particular form to fit the data). 
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quoted form for this relation, ignoring the weak mass de- 
pendence. Remarkably, our general conclusions should hold 
for haloes in any mass range 



2.2 Generalisation of the standard formulation 

The standard approach to relaxation in gravitational sys- 
tems has developed along the foundation set by Chan- 
drasekhar (1943). In analogy with the case of dilute gases, 
it is assumed that perturbations experienced by a test par- 
ticle arise from independent and local two body encounters. 
Since the mean field along a particle's trajectory changes on 
a timescale comparable to its dynamical time td, while the 
timescale of an average encounter is ~ td/N, this assump- 
tion is generally satisfied; encounters with duration ~ td 
being relatively rare. As the number of particles in the rep- 
resentation of a self-gravitating system is increased, the ef- 
fect of strong encounters becomes increasingly unimportant; 
in softened systems, the subject of cosmological simulations, 
strong encounters are completely suppressed if the softening 
length is of the order of the interparticle distance. The stan- 
dard formulation of relaxation theory therefore assumes that 
the effect of discreteness noise can be modelled in the form 
of weak random encounters. Numerical tests (e.g., Huang, 
Dubinsky & Carlberg 1993), seem to vindicate these basic 
premises. 2 

The product of weak, local and independent encoun- 
ters is a diffusion process, whereby a particle's dynamical 
variables undergo a random walk around their unperturbed 
values. For a spherical system with isotropic velocities that 

2 Note nevertheless that such work generally focused on energy 
changes along particle trajectories. The response of the trajec- 
tories themselves (and thus the effect on quantities like angular 
momentum) can exhibit significantly stronger sensitivity to noise 
— even in the simplest case of deflections caused by fixed back- 
ground particles (Athanassoula et al. 2001; El-Zant 2002). 



can approximated in terms of a Maxwellian, the mean dif- 
fusion coefficient within radius r 

2 = G^nA 
Ka 

describes the average rate at which the square velocities of 
particles deviate from those of unperturbed trajectories in 
a corresponding smooth system (with iV — > oo). Here A 
is the ratio of maximum to minimum impact parameters 
and K = 1/15.4 (according to Spitzer & Hart 1971). In 
quasiequilibrium, at any given instant, the average velocity 
dispersion inside radius r can be calculated from 

(j 2 (r) = — / Vrdr = 7 2 u 2 = j 2 GM/r = 7 2 G —irpr 2 . (9) 
r Jo 3 

The local one dimensional velocity dispersion i> 2 can be eval- 
uated by solving the Jeans equation (cf. Appendix B). The 
above relations thus define the weak function of radius 7 that 
is used in the calculations below. This definition should hold 
because the virialised region of a cosmological halo remains 
near equilibrium through most of its evolution — major 
mergers being relatively rare, and even when they do oc- 
cur the system is out of equilibrium for a time td generally 
much smaller than the (relaxation) timescales of interest. 

When the system parameters are changing, as during 
the growth of a cosmological halo, the mass inside any given 
radius will be time dependent. As a consequence, there will 
be corresponding variations in p — ^3^3 and 7. Neverthe- 
less, if the timescale between encounters giving rise to the 
relaxation process (r e ~ td/N) is much smaller than the 
timescale for variation of the system parameters r s , then 
local (in time) averages are allowed, and the diffusion coeffi- 
cient becomes a well defined function of time. From Eq. JSJ, 
and because 

\z\ = (l + z)H (10) 
and 
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pc=V 



3H 2 
8^G' 



(11) 



one finds that t s ~ -rr— 

* 1 + 2 



S7r Q p ■ Therefore, unless z is very 
1/^/vGpc (recall that v > 200); with the 



large, t s > td 

implication that the temporal locality condition (r e <C r s ) 
is in fact weaker than that for spatial locality (r e <C ru), 
required for the validation of the diffusion approach. 

The expression for the mean diffusion coefficient, inside 
radius r and at time t, can then be written as 



{(Av) 2 )(r,t) _ y/Gp(r,t) 



a 2 (r,t) Kj 3 (r,t)-y/4ir/3 M(r,t) 



lnA(r,t). (12) 



The relative mean square perturbation to particle velocities 
due to encounters with other particles during the halo evo- 
lution is 



(vl/a 2 ) = 



((AiQ 2 )M) 

a 2 (r,t) 



<it. 



(13) 



where tf refers to some chosen initial, reference, formation 
time (for any given r measured at z = 0, we will take the 
redshift to be the time when r v - u (tf) = r) and to refers to 
the end of the simulation (assumed to correspond to z — 0). 

Using O and JT^J one can define the relaxation time 
in a simulated cosmological halo implicitly: 



/Up 



M 



In Adt 



(14) 



If treiax < to the dynamics can be expected to be completely 
dominated by discreteness noise. This is in line with the 
standard application to time independent systems, where 
tf = and {v 2 /a 2 } is assumed to be time independent; the 
the result are the familiar expressions for the relation time 



tr. 



K 



ty M 



G 2 pm In A 



lnA^/Gp m 



N 

°- 1 InA TD 



(15) 



(where N(r) and td{t) are the particle numbers and mean 
dynamical time inside r). 

To fix the value of A, we assume that the resolution 
of the simulations in question corresponds to the local in- 
terparticle distance. The average of this quantity within ra- 
dius r will define the minimum impact parameter 6 m i n = 
(|'7rr 3 /A r (r)) 1/ ' 3 . The maximum impact parameter will be 
taken to correspond to the virial radius. We therefore have 



A(r,t) = ^= (±) 1/3 N^(r,t) 

Omin V 4-7T / 



(16) 



3 SIMPLE ESTIMATES 

In the next section we will numerically evaluate the integral 
in Eq. 1141 . invoking the time dependence appropriate in 
the currently favoured ACDM cosmology. Some insight can 
however be gained by defining characteristic relaxation times 
for an evolving cosmological halo and for its substructure. 
This is the subject of this section. 



3.1 A characteristic relaxation time for 
cosmological haloes 

We consider first the evolution of the parent halo, or most 
massive progenitor, and assume that all mass is added to it 
in the form of a smooth component. Of course, in reality, 
accretion of material takes the form of clumpy subhaloes, 
but most of their mass is rapidly stripped, so that, at any 
given time, there is usually only 10 — 20% of the mass of the 
halo in the substructure. Their internal relaxation is dealt 
with separately below. 

If the relaxation time within some radius r is smaller 
than the timescale r s for the mass distribution to signifi- 
cantly evolve, it is likely that the system is heavilly affected 
by discreteness noise within radius r — the rationale being 
that if the local relaxation time is small compared to the 
time it takes for more particles to be added to the system 
(thus decreasing the relaxation rate) then discreteness noise 
will have a significant effect. 

A local relaxation time can be obtained by freezing the 
system at some time t = t(z) and using (1151 . Comparing the 
two timescales one can write 



R — ^rclax / 7"s 



iferyVt 71 " \dz/dt\ M 



InA v/Gp 

and eliminate dz/dt by invoking fj). (1101 and (II I li to get 

0.087 x /27^7rQ7 : 



(17) 



InA 



-(l+z)N(r, 



(18) 



where (by Eq.|5J the number of particles inside r at redshift 
z can be expressed in terms of the final number of particles 
with which the halo is resolved A^o = Mo/m: 



No e~ 



(19) 



The mass ratio entering into this last relation can be calcu- 
lated via either Eq. (0 or Eq. 

The dimensionless relaxation time R has to be signifi- 
cantly greater than unity, at all z, if the simulated halo can 
be considered free of artificial relaxation inside radius r. 3 
What is required therefore is that 



A?o > 11.5 



InA 



Mo 
M 



(20) 



Tray ivi 1 + z 

In Fig. [5] we show the variation of the of the quantity 
(t[tx)No as a function of radius at different redshifts. The 
radii are expressed in terms of the final virial radius (at 

2 = 0), by transforming them using Eq. EJ. The cutoffs m 
the curves correspond to the virial radius at the denoted 
redshift (that is the maximum size of the halo considered at 
that redshift). 

From this figure it is apparent that, provided 7 3 ~ InA, 

3 Note that, because the virial radius and concentration are time 
dependent, accreted mass is not deposited uniformly; there is a 
preference, especially at later times, for mass increase in the outer 
regions (e.g., Zhao et al. 2003). However this implies that the char- 
acteristic time for mass change in the inner regions, where most 
of the relaxation occurs, is larger than that predicted by Eq. ISt . 
The condition _R > 1 thus constitutes a minimal requirement. 
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Figure 2. Number of particles No required to compose a final main halo at z = so that, at indicated redshifts, the mean relaxation 
time will be larger than the characteristic time for the the halo to increase its mass at that redshift (see Eq. |§J. Radii are rescaled to 
units of the virial radius at z = 0. The plot on the left corresponds to NFW haloes while that on the right to the M99 profile. Note that 
7 3 /lnA ~ 1 



that artificial relaxation may have a dominant effect inside 
the inner ~ 1% of the final virial radius for N ~ 10 6 . Note 
that the demand for larger particle number is most stringent 
at larger times (smaller z). This is because the characteristic 
growth time r s is a steeply increasing function of z. The 
figure also suggests that, for fixed 7, the relaxation effects 
are more pronounced in the inner regions of NFW haloes, 
compared to the M99 case. Nevertheless, as we will see in 
Section [4.21 the introduction of variation in 7 reverses this 
situation, since 7 is significantly larger in the inner regions 
of NFW systems (this follows from the variations of the ratio 
of velocity dispersion to circular speed shown in Fig. Bl). 




3.2 Characteristic relaxation time for 
substructure 

We consider a subhalo accreted at redshift z a and remain- 
ing a separate dynamical system up to z = — i.e., it stops 
growing, its outer regions are in fact stripped, but it keeps 
a dynamically distinct core. For this purpose we exploit the 
possibility of relating the relaxation and the local Hubble 
time; and in order to get simple closed form estimates, we 
will neglect relaxation at times t < t(z a ) and focus on relax- 
ation effects since accretion. From equations 1151 . 1101 and 
lllll we can define a characteristic relaxation time 



&rela 



0.0877r 7 d 



(21) 



where M, v, 7 and A are all measured inside radius r at 
z = z a - This timescale should be compared to the time in- 
terval separating the accretion epoch t(z a ) and the end of 
the simulation at to = t(0). This determines a characteristic 
substructure relaxation time given by 



tt. 



c(r, Z a ) 



to — t{Za) 



(22) 



10" 



10 



Figure 3. Relaxation time of substructure haloes consisting of 
a thousand particles accreted at redshift z a . The timescales are 
expressed in terms of the time tjj elapsed between epochs z a 
and z = and are given at different fractions of the final virial 
radius. Values smaller than unity correspond to situations where 
the dynamics is completely dominated by artificial relaxation. 
Haloes are assumed to be of the NFW type 



Now further assume that t(z) — |_ff and z = 
(to/t(z)) 2 ^ 3 — 1, as is appropriate for an Einstein de-Sitter 
universe (in the next section we verify our results by under- 
taking the full calculation, including pre-accretion evolution, 
in the currently favoured "concordance" cosmology 4 ). In this 
case we have 



4 Note that while t = t(z) for ACDM is quite different from that 
in an Einstein-de-sitter model the ratios t(zi)/t(z2), entering into 
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0.087tt7 3 



N(r, 



In A 



(Za + 1) 3/2 



(23) 



If, within some given radius, this quantity is less than unity 
the implication is that, by z = 0, the dynamics has become 
dominated by artificial discreteness noise. 

In Fig. [3] we reproduce several plots where this dimen- 
sionless relaxation time, thus defined, is shown for haloes 
assumed to have N(r v i T ,z a ) = 1000 when they are sub- 
sumed into their parent structures at different z a (in the 
case of satellites haloes we only reproduce the results for 
NFW haloes, those for the M99 haloes are very similar). 
As can be seen, especially at relatively small radii, the re- 
laxation time can be significantly smaller than the subhalo 
lifetime. Furthermore, it is to be noted that it will be this 
inner core of the subhalo that will survive stripping (see also 
Fig. H]] for more detail concerning the radial distribution of 
the relaxation time and Section f5.2l for a discussion of some 
possible consequences). 



4 DIRECT CALCULATIONS 
4.1 General considerations 

From equations 11 21 and I|1H[I . the expected relative mean 
square perturbation due to discreetness imposed on to a test 
particle moving within an evolving cosmological halo is 



(v p /a ) 



2K 



taA r,t) v ,, / \ T , K ! dt. (24) 
v ' ; 7 3 (r,t) N(r,t) v ; 



The integral is evaluated within a given fixed fraction of 
the final virial radius r vir (to) — with the consequence that 
the radial coordinate contains an implicit time dependence 
(i.e, r — r(t)). The relations in Section [2.11 can be used to 
determine this dependence, as well as that of the other quan- 
tities entering into 12411 . provided the evolution of redshift 
as a function of time is known. For a flat universe with cos- 
mological constant, the required transformations are given 
in Appendix A (in the calculations below we use Qy = 0.7 
and h — 0.7). Since the virial radius increases with time, a 
given fraction of the virial radius at t — to will correspond to 
the whole virial radius at some earlier time. This naturally 
determines the lower limit of integration, the formation time 
tf — for some fixed fraction of r v i r (z — 0), it corresponds 
to a redshift zj where this fraction is equal to whole virial 
radius r v it{zf). 



4.2 Relaxation of the most massive progenitors 

In Fig. 0] we show the RMS perturbation due to particle 
noise, as a function of radius, expected for a halo with mass 
evolving according to Eq. (with a = 1) for several final 
values of the final total particle number No. They are ob- 
tained by solving 1241 using an adaptive integrator (NAG 
D01AJF) with a tolerance of 10~ 4 . 



the equation below, differ by a factor of at most by ~ 3% for 
10 > z > 1 and by ~ 30% up to z = 0. 



An interesting characteristic of the plots in Fig. |l]is the 
perfect power law behaviour of (Vp/a 2 ) 1 ^ 2 over most of the 
radial interval. This is a consequence of a curious property 
of cosmological haloes, already noted by Taylor & Navarro 
2001; namely the power law form of the phase space density. 
It is this phase space density, ~ p/o" 3 , that determines the 
rate of relaxation. This property enables one to deduce a 
particularity simple fit for the variation of the relative RMS 
perturbation: 



2\l/2 



10" 




(25) 



where r/r v i r refers to the fraction of the virial radius at 
z — 0. The number of relaxation times inside radius r can 
be counted as follows: 



Tl (trelax) = ("p/o" 2 } ~ 10 



10 

N 



V T ) ' 



(26) 



which of course needs to be much smaller than one if artifi- 
cial relaxation is to be negligible within radius r. 

For No < 10 6 then, particle motion is expected to be en- 
tirely dominated by noise in the inner percent of the virial ra- 
dius; regions bounded by radii an order of magnitude smaller 
still having undergone ~ 100 relaxation times. These results 
are in agreement with the simple estimates presented in Sec- 
tion ljMl fwith the difference that the effect on M99 haloes are 
larger here, because variations in 7 is taken into account). 

The inner 0.1% is one relaxation time old even for No = 
10 s . Indeed, a final particle number inside the halo's virial 
radius of No = 10 12 is required for noise to be reduced to 
the level of a few percent there. For a million particles this 
noise level is only reached at almost a third of the final virial 
radius; and this is larger than the average halo scale length 
r s = r vil /c (c.f. Eq.0. 

Fig. shows, for the inner percent of the final virial 
radius, the z-evolution of the relaxation rate (the integrand 
in 1241 along with the principal components determining it. 
It shows that the relaxation rate decreases with redshift, for 
both the NFW and M99 profiles — the effect being more 
pronounced in the former case because the ratio of velocity 
dispersion to circular velocities that determines 7 decreases 
faster with radius (Appendix B); and at smaller z because, 
for a given fixed fraction of the final virial radius, smaller 
radii are probed (recall that at Zf the whole virial radius 
corresponds to the inner one percent at z = 0). 

Finally, note that, since new particles are continuously 
being introduced into the region, the way in which particles 
are affected by relaxation will differ, with the probable con- 
sequence that relaxation driven evolution can be expected 
to be different than in a static NFW type system. 



4.3 Relaxation of subhaloes 

In Fig. |S| we show the RMS perturbation in velocities ex- 
pected as a result of particle noise for a subhalo composed 
of a thousand particles when it stops growing (after being in- 
corporated into a larger structure). Before this, we assume 
that its mass grows at the exponential rate described by 
Eq. J5J with (as always) a = 1. These results are in close 
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Figure 4. Expected RMS perturbation to the velocity resulting from particle noise in evolving most massive progenitor cosmological 
haloes. It is evaluated at a given fraction of the final virial radius by integrating Eq. 1241 . from the time when that fraction was equal 
to the whole virial radius at epoch t = tf up to z = (t = to). Results for NFW haloes are shown on the left while those corresponding 
M99 ones are on the right 



0.1 x N (0.01 r vjr (z = 0) ) 
- i n m i 1 1 i h i ii i n i i n 





Figure 5. Instantaneous values of the perturbation to the velocities that arises from discreteness noise inside 1% of the final virial radius, 
and the corresponding magnitude of terms figuring in Eq. (eq:msdirect) that contribute to this quantity for NFW haloes (left) and M99 
profiles. All contributions are of similar magnitude for the two profiles, except for the ratio of velocity dispersion to circular velocity, 
which is larger at small radii (hence lower redshift) in the NFW case. This leads to a somewhat smaller perturbation 



agreement with the simple estimates derived in Section f3. 21 
they reinforce the conclusion that relaxation can be a signif- 
icant factor determining the structure of subhaloes in simu- 
lations, their distribution and their effect on the parent halo 
structure (Section l5.2l l. 



5 CONCLUSION 

This paper presented an attempt at assessing the impor- 
tance of particle noise in the evolution of collapsed cold dark 
matter structures, by generalising the diffusion formulation 
first proposed by Chandrasekhar to the case when the mass 



distribution and radial extent of the system under consider- 
ation are variable. Although a most massive progenitor of a 
z = halo will typically have increased its collapse mass by 
more than an an order of magnitude by accreting subhaloes, 
at any given instant > 80% of its mass is in a smooth com- 
ponent, most material in subhaloes having been stripped 
(e.g., Gao et al. 2004). This enables one to divide the anal- 
ysis in two parts. In this context, 'parent' haloes, i.e ones 
characterised by continuously increasing (virial) mass and 
radius throughout their evolution, were grown according to 
the empirical formula discovered by Wechsler et al. (2002). 
Subhaloes were treated differently, by assuming that their 
growth is arrested at accretion. Stripping was not explicitly 
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Figure 6. RMS perturbations to the velocity of satellite haloes 
composed of a thousand particles when they stop growing due to 
incorporation into a larger structure at the redshifts shown (up 
to which they are assumed to grow in mass, and hence particle 
number, according to Eq.|5J. Results are obtained using Eq. 1241 
and correspond to subhaloes with NF W profiles. Horizontal lines 
show perturbations corresponding to (from bottom up) one and 
ten relaxation times 

taken into account. But since it is in the inner region of sub- 
haloes that relaxation is most significant, this does not ap- 
pear to represent a major idealisation — as high resolution 
simulations, specially designed for the purpose of examining 
the issue, show the inner regions of stripped haloes to be 
largely unaffected (Kazantzidis et al. 2004). 

We have derived simple closed form formulas for the 
characteristic relaxation times (Section [3J, and also inte- 
grated the relevant equation characterising the expected 
RMS perturbation to the velocity directly (Section |1J, The 
predictions of the two approaches agree quite closely. In 
what follows we summarise our results and sketch some pos- 
sible consequences. 

5.1 Relaxation in parent haloes 

Our results were not found to be significantly affected by the 
exact form chosen for the evolution of the concentration pa- 
rameter c or the parameter a characterising the mass growth 
rate (cf. Eq. both quantities that are dependent on the 
final halo mass (in physical units), so this latter quantity 
does not enter directly into our presentation here, which is 
concerned only with the relaxation of particles after they 
have been accreted onto the most massive progenitor and 
are subsequently part of its smooth component (in the next 
two subsection we consider relaxation inside substructure, 
and briefly comment on the expected mass dependence of 
noise-induced relaxation in terms of merger history). 

The power law (with index ~ — 2) radial variation of the 
phase space density of cosmological haloes ensures a steep 
dependence of the relaxation rate, a linear function of this 
quantity, on radius. The perturbation due to discreetness 



noise is adequately quantified by the empirical fits given in 
Section ET31 

A test particle that is present inside the inner percent of 
the final virial radius, from the point in time when the virial 
radius of the halo was equal to this fraction, will experience 
a relative RMS perturbation to its velocity ~ 1 when the 
final halo is resolved with 10 6 particles — implying that its 
dynamics is completely corrupted by discreteness noise. In 
the very inner few thousandth of the final radius the RMS 
perturbation due to discreteness noise is an order of magni- 
tude larger, meaning that these regions are of the order of 
a hundred relaxation times old. 

In an evolving cosmological halo, particles are continu- 
ally added within any given radius; it is also likely that some 
can gain energy (e.g. by interaction with sinking satellites; 
cf. El-Zant et al. 2004, El-Zant 2005) and move into trajec- 
tories with larger average radii. Particles that are accreted 
at a later stage are more mildly affected by relaxation. At 
any given fraction of the final radius therefore, there will be 
a range in the severity of artificial relaxation, depending on 
the accretion epoch of individual particles. This somewhat 
complex situation implies that effects of relaxation may be 
rather different than the familiar situation of fixed-mass sys- 
tems, and therefore difficult to detect. 

Our calculations thus represent a quantification of the 
magnitude of the perturbation of the discreteness noise, 
rather than its effect on the detailed dynamics and macro- 
scopic structure. While it is quite probable that the signifi- 
cant role for artificial relaxation predicted by our analysis is 
reflected in the contradictory claims concerning the conver- 
gence and shape of the inner density profile of CDM haloes 
(Navarro et al. 2004; Diemand et al. 2005), with the pre- 
dicted effect for static systems being a weakening of the 
cusp, further work is needed to determine the consequence 
in the case of evolving cosmological haloes. 

The RMS perturbation to the velocities decreases as 
~ 1/yN, implying that a parent halo needs to be resolved 
with N > 10 10 particles at z = if discreteness noise is not to 
dominate the dynamics of the inner 0.1% of the virial radius; 
and if it is to be negligible ((Vp/a 2 ) 1 ^ 2 < 1%) in the whole 
region where r > 0.1r a ~ 0.01r v i r - While perturbations of or- 
der unity are likely required for significant energy relaxation 
and restructuring of the azimuthally averaged density pro- 
file, much smaller perturbations may affect particle trajec- 
tories, with accompanying modification (isotropisation) of 
the velocity distribution and loss of triaxiality in spatially 
asymmetric systems (e.g., Merritt & vallurri 1996; see also 
Section l2~2jl . 5 



5 The reason this happens is due to a transition from regular 
to chaotic motion, which is very sensitive to perturbation. A very 
simple example of this phenomenon is the case of a pendulum on a 
'separatrix' trajectory passing near the unstable (upper) equilib- 
rium. Small perturbations can turn (regular) trajectories rotating 
in one direction into ones rotating into the opposite direction, or 
oscillating without a definite sense of rotation, or even transit 
between all these different possibilities (chaotic trajectories). 
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5.2 The relaxation of substructure 

The results presented in this paper suggest that the situation 
with substructure is still more severe. A parent halo that is 
identified with No ~ 10 7 has (from Eq. ~ 10 6 particles 
at z ~ 2. Suppose a subhalo with N su b = 0.001AT is then 
accreted. Assume that between redshifts two and one this 
1000-particle halo will be stripped of ~ 90% of its mass (as 
would be expected from Gao et al. 2004 Fig. 14). For an 
NFW halo this also corresponds to a truncation of about 
90% in radius. But according to Fig. [3] this inner region will 
be significantly affected by relaxation. In fact, by z — 0, it 
is expected to have undergone enough relaxation times to 
have approached core collapse! But this fate is only likely if 
it is not, in the meantime, dissolved by stripping. 6 

What effect is the significant relaxation expected to 
have on stripping? Cosmological haloes have a 'tempera- 
ture inversion' in their inner velocity distribution (cf. Ap- 
pendix B). Relaxation therefore causes initial expansion of 
the inner region, forming an isothermal core, before it re- 
contracts like a relaxing globular cluster (e.g., Hayashi et al. 
2003). In the first phase, the core becomes less dense and 
more easily stripped; in the second the situation is reversed. 
Stripping is most efficient for subhaloes that venture near 
the centre of the parent, experiencing strong tidal forces. 
If, while the outer regions of a subhalo are stripped dur- 
ing successive passages, the inner region relaxes to a more 
diffuse density distribution, further stripping will be accel- 
erated and the core may completely dissolve. Conversely, if 
the stripping is slow, there may be sufficient time for core 
contraction to take place. Further stripping is subdued and 
core dissolution becomes less likely. 

If, as is suggested in several studies (e.g., Syer & White 
1998; Dekel et al. 2003b; Dekel et al. 2003a), the structure 
of the halo profiles found in numerical simulations is de- 
pendent on the interaction between the subhaloes and their 
parents, then it is clear that such significant re-engineering 
of the subhaloes can be of major importance in determining 
the parent halo profile. It obviously also has implications for 
the number and spatial and orbital distribution of substruc- 
ture. Given the discussion above, one would expect artifi- 
cial relaxation to enhance the number of haloes that do not 
venture into the inner regions (i.e., those accreted on less 
eccentric orbits) and decrease the number of those that do 
venture there. 



which have spent considerable time inside subhaloes before 
being stripped. 

Simulations suggest that substructure has a mass func- 
tion such that dn/dm ~ m~ 9//5 , quite independently of red- 
shift (e.g., Gao et al 2004). If one ignores the logarithmic 
dependence, the relaxation time when expressed in units of 
dynamical time is, for any given subhalo, a linear function 
of the number of particles in it. Consider then a mean re- 
laxation time averaged over subhaloes and normalised over 
an averaged dynamical time, 



(tra) 



N f 

J rr. 



- g / 5 dm 



f 

■J rr. 



m~ g / 5 dm 



(27) 



where the integration limits correspond to the least and most 
massive subhalo present. When the latter is orders of mag- 
nitudes more massive than the former, the above relation 



4/8 



If both and m max are in- 



gives (t ra ) ~ N m max m. 
dependent of iVo then the relaxation time scales linearly in 
with No as in the standard case. 

If, on the other hand, one supposes that, because of in- 
creasing resolution, m m i n ~ 1/No (proportionately smaller 
haloes appear in the resolved field as No increases) then 
(t ra ) ~ Nq , which is roughly the scaling claimed by Die- 
mand et al. (2004) on the basis of direct calculations of the 
relaxation time in the context of cosmological simulations. 
Their results also show that the N-scaling is closer to the 
canonical linear relation as one moves towards the centre of 
the main halo. This is also to be expected; since, near the 
centre, stripping causes the fraction of mass in subhaloes to 
decrease, while the relaxation in the smooth (parent halo) 
component becomes more efficient; and this scales roughly 
linearly with particle number. 

Finally, note that although we have not followed the 
detailed merger history, the expectation derived from the 
work presented here is that relaxation is more enhanced in 
cluster haloes — since these form relatively recently, from 
particles that spend most of their existence inside poorly re- 
solved structures. In contrast, a galaxy sized halo acquires 
most of its mass (i.e., particles) at larger redshift. It fol- 
lows that, given the same resolution at z — 0, cluster halo 
particles would have, on average, been more affected by dis- 
creteness noise. This is in line with the conclusion reached 
by Diemand et al. (2004). 



5.3 The N-scaling of the substructure-dominated 
relaxation time 

The expectation from the results just outlined is that, be- 
yond the inner percent or so of the final virial radius, relax- 
ation will be dominated by particles inside substructure, or 
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6 Note that, because relaxation is quite a sensitive function of 
radius (i r dax ~ r2 from Fig. [S] and Eq. 1261 this same argument 
can be made for more central regions of more massive subhaloes 
(e.g., for the inner 1% of a subhalo consisting of 10% the mass of 
mass of the parent halo, accreted at z ~ 2 when the parent halo 
had N ~ 10 6 ), or more extended regions of less massive ones. 
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Figure Bl. Velocity dispersion in spherical cosmological haloes 
with isotropic velocities in terms of the local circular velocities as 
a function of radius scaled to the halos scale length r s . 
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In a similar manner we find the critical density to vary as 
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Using these formulas the integral in Eq. 124H can be evalu- 
ated with respect to the variable Cvt. When the constant in 
Eq. 1A4I is evaluated by assuming z = at t = to, the result 
does not depend on the value of to- This is analogous to the 
familiar case whereas two body relaxation in an isolated sys- 
tem does not depend on the time elapsed in physical units, 
but instead on the number of dynamical times the system 
has been through. 



APPENDIX A: TIME DEPENDENCE OF p c 
AND Z IN FLAT UNIVERSES WITH 
COSMOLOGICAL CONSTANT 

According to Kolb & Turner (P. 55), for a flat universe with 
cosmological constant one has 



t = -h~ l q 
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which, using fV = 4h and C v 
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which by integration (with A the associated constant) gives 



APPENDIX B: VELOCITY DISPERSIONS 

The Jeans equation for spherical system with isotropic ve- 
locity dispersion can be written as (e.g., Binney & Tremaine 
1987) 



d(pVr 

dr 



dr ' 



(Bl) 



where v% = |a 2 is the radial velocity dispersion. The general 
solution is 



(B2) 



For any p — > 0, as r — > oo, one must have C = 0, if the 
velocity dispersion is to be bound at large radii. The ratio 
of this unique4 physical solution to the local circular velocity 



(B3) 
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the form reflecting the fact that, as opposed to circular mo- 
tion, the velocity dispersion at any r is caused by particles 
that venture in and out of that radius. 

For power law density distributions p a r"™ solutions 
of Eq. lIHlJ, with C = 0, imply that 

„72 1 

(B4) 



v 2 c 2n - 2 ' 
for n/1 and 

4= -Inr, (B5) 

Vc 

(r < 1). Note that, in the central regions, vl oc — rlnr for 
the NFW profile and oc r 1/2 for the M99 profile. Spherical 
cosmological haloes thus, necessarily, have a 'temperature 
inversion', at least as long as the velocity dispersion can be 
considered isotropic. 

In Fig. IB1I we reproduce the ratio of the three dimen- 
sional velocity dispersion to circular velocity for both profiles 
dealt with in this paper. It determines the value of gamma 
(corresponding to the average of this quantity within any 
given radius) that enter into the equations calculating the 
effects of relaxation. 
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